Constrained CPD of Complex-Valued Multi-Subject fMRI Data via Alternating Rank-R and Rank-1 Least Squares

Complex-valued shift-invariant canonical polyadic decomposition (CPD) under a spatial phase sparsity constraint (pcsCPD) shows excellent separation performance when applied to band-pass filtered complex-valued multi-subject fMRI data. However, some useful information may also be eliminated when using a band-pass filter to suppress unwanted noise. As such, we propose an alternating rank-R and rank-1 least squares optimization to relax the CPD model. Based upon this optimization method, we present a novel constrained CPD algorithm with temporal shift-invariance and spatial sparsity and orthonormality constraints. More specifically, four steps are conducted until convergence for each iteration of the proposed algorithm: 1) use rank-R least-squares fit under spatial phase sparsity constraint to update shared spatial maps after phase de-ambiguity; 2) use orthonormality constraint to minimize the cross-talk between shared spatial maps; 3) update the aggregating mixing matrix using rank-R least-squares fit; 4) utilize shift-invariant rank-1 least-squares on a series of rank-1 matrices reconstructed by each column of the aggregating mixing matrix to update shared time courses, and subject-specific time delays and intensities. The experimental results of simulated and actual complex-valued fMRI data show that the proposed algorithm improves the estimates for task-related sensorimotor and auditory networks, compared to pcsCPD and tensorial spatial ICA. The proposed alternating rank-R and rank-1 least squares optimization is also flexible to improve CPD-related algorithm using alternating least squares.


1) Simulated fMRI data analysis
We examine noise effects of different SNR levels on ARRR1LS with several CPD algorithms (i.e., accALS, enALS, AO-ADMM, NLS, MINF, and COMFAC) under simulated fMRI data with spatial change, as shown in Fig. S1. Note, since these CPD methods do not incorporate the shift-invariance property, we here do not consider the temporal time delay differences. We can conclude from Fig. S1

2) Actual fMRI data analysis
In order to further show the outperformed separation performance of ARRR1LS, we present

C. Results of Task-related Auditory Component
For the experimental fMRI data, another task-related component is auditory component. Thus we further present the comparison results of sARRR1LS-PO, T-sICA, and pcsCPD in this subsection. Since group ICA is widely-acceptable for multi-subject fMRI data analysis and magnitude-only fMRI data persist smaller noise than complex-valued fMRI data, we use auditory SM extracted by group ICA [S13] of magnitude-only multi-subject fMRI data as the auditory SM reference in Fig. S2(a). Besides, we generate the task-related auditory TC references by convolving the stimuli with the canonical SPM hemodynamic response functions in Fig. S2(b).  Fig. S3, pcsCPD is higher than T-sICA in filtered fMRI data analysis but is lower than T-sICA in raw fMRI data analysis. This indicates that pcsCPD is sensitive to noise when extracting auditory component. Moreover, as the auditory component is more difficult to decompose than sensorimotor component, these three methods all show lower average  values of auditory shared SM and TC estimates in noisier raw fMRI data analysis than in filtered fMRI data analysis. In addition, these three methods comprehensively obtain the larger standard deviations of  values in raw fMRI data analysis than in filtered fMRI data analysis as shown in Fig. S3. We also present Fig. S4 to show magnitude and phase parts of the typical auditory shared SMs estimated by T-sICA, ICA-sCPD, pcsCPD, and sARRR1LS-PO in raw and filtered fMRI data analyses. The Vall, Vin, Vout, and Vin/Vall values of typical sARRR1LS-PO has the smallest noise voxels (i.e., smallest Vout  Table SII. From Fig. S4 and Table SII, the shared SMs estimated by proposed and right middle, superior and inferior temporal areas middle, superior and inferior temporal areas in both two analyses. We can conclude that these three methods relatively obtain better shared SMs (i.e., higher  values and Vin values) in filtered fMRI data analysis than in raw fMRI data analysis. The pcsCPD extracts better shared SMs (higher  , Vin, and Vin/Vall values and lower Vout value) than T-sICA in filtered fMRI data analysis, while obtains the worst shared SMs in nosier raw fMRI data analysis (see Fig. S4 and Table SII). For auditory component, we give Fig. S5 to present the results of shared TC magnitude and phase parts, subject-specific time delays and intensities in raw and filtered fMRI data analyses. The magnitude and phase parts of shared TCs estimated by the proposed sARRR1LS-PO can better follow fluctuations of auditory TC reference than other two methods in both two cases in Figs. S5(1)~(2). Moreover, the number of correctly estimated time delays of the proposed sARRR1LS-PO is larger than pcsCPD as shown in Fig. S5(3). More specifically, the sARRR1LS-PO not only can accurately estimate the larger time delay (e.g., subject 1), but also can accurately estimate the zero time delays. Furthermore, subject-specific intensities curves of pcsCPD are closer to those of sARRR1LS-PO than those of T-sICA in both two analyses (see Fig. S5(4)). More specifically, the  values are 0.824 vs. 0.718 in raw fMRI data analysis and 0.924 vs. 0.853 in filtered fMRI data analysis.

D. Comparison Results of Orthonormality Constraint
In this subsection, in order to evaluate the importance of orthonormality of the proposed method, we further evaluate the separation performance of the proposed method without orthonormality constraint (shorted as sARRR1LS-P) and with orthonormality constraint (i.e., sARRR1LS-PO) as show in Table SIII. As for the raw fMRI data analysis, sARRR1LS-PO shows obviously higher means and lower standard deviations of ρ values of sensorimotor shared SM magnitude, SM phase, TC magnitude, and TC phase estimates than sARRR1LS-P. In the filtered fMRI data analysis, sARRR1LS-PO also shows higher average ρ values than sARRR1LS-P. However, the ρ value differences between sARRR1LS-PO and sARRR1LS-P in filtered fMRI data analysis are smaller than in raw fMRI data analysis. This indicates that the orthonormality constraint can effectively improve the performance of the proposed method for noisy raw complex-valued fMRI data analysis.

E. Comparison Results of Computation Cost
We finally analyze the computation cost of the proposed sARRR1LS-PO method with other methods from the following aspects: computation complexity, computation time and maximum computation memory for actual fMRI data experiment. First, we compare the computation complexities of updating shared SMs S, shared TCs B and subject-specific intensities C for the proposed method and pcsCPD, as shown in Table SIV. The complexity of the proposed method is approximately R/V 3 of that of pcsCPD. Since R<<V (i.e., R = 50 vs. V = 59610 in actual fMRI experiment), the proposed method has much lower computation complexity than pcsCPD.  Second, we illustrate the computation time and maximum computation memory for T-sICA, pcsCPD, and sARRR1LS-PO in Table SV. Notice that we run the methods on MATLAB2020a under Win10 operating system with Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz and 32G memory.
Although T-sICA costs the minimum computation time and memory, it shows the worst separation performance for both simulated and experimental fMRI data analyses (see Figs. 3 and 5), compared with pcsCPD and sARRR1LS-PO. The proposed sARRR1LS-PO requires 702.2% lower computation time and 178.3% lower computation memory than pcsCPD. Above all, the proposed method not only shows better separation performance but also owns lower computation complexity and higher speed than pcsCPD.